

library(tseries)
library(MASS)
library(forecast)

myplot = function(X, title="", lag.max=NULL) {
  layout(matrix(1:3, 3, 1))
  plot.ts(X, main=title)
  acf(X, main=title, lag.max=lag.max)
  pacf(X, main=title, lag.max=lag.max)
}

## study of seasonality of deaths

deaths
myplot(deaths)

# the 12-differencied series
deaths.diff = diff(deaths, 12)
myplot(deaths.diff)

# any other lag shows seasonality
showlags = function(X) {
 layout(matrix(1:16, 4, 4, byrow=TRUE))
 for(d in 1:16) {
   acf( diff(X, d), main=paste("lag = ", d))
 }
 quartz()
 layout(matrix(1:16, 4, 4, byrow=TRUE))
 for(d in 1:16) {
   plot.ts( diff(X, d), main=paste("lag = ", d))
 }
}
showlags(deaths)

## spectrum inference for deaths
layout(matrix(1:2, 2, 1))
X <- deaths
spectrum(X) # the periodogram
m <- 3
spectrum(X, spans=c(m, m)) # the smoothed periodogram
# on trouve un pic la ou il faut...

X <- deaths.diff
layout(matrix(1:2, 2, 1))
spectrum(X) # the periodogram
m <- 3
spectrum(X, spans=c(m, m)) # the smoothed periodogram

# band = (2*m+1)*(2*pi)/length(deaths)*1/sqrt(12)

### a fit for lh data

X <- lh
myplot(X)

n <- 1000
ar1 <- arima.sim(n=n, model=list( ar=c(0.6, -0.3 )))
myplot(ar1)

# on fitte un AR(1)
lh.ar1 <- arima(X, order=c(1, 0, 0))
lh.ar3 <- arima(X, order=c(3, 0, 0))
# ar(lh, aic=FALSE, order.max=1)
tsdiag(lh.ar3, gof.lag=30)
# find the best fit according to AIC
lh.bestfit <- auto.arima(X)
quartz()
tsdiag(lh.bestfit)
# prediction
lh.fore <- predict(lh.bestfit, 12)

# intervalle de prédiction
ts.plot(lh, lh.fore$pred, lh.fore$pred + 1.96*lh.fore$se, 
        lh.fore$pred - 1.96*lh.fore$se, 
        gpars = list( lty=c(1, 1, 2, 2), lwd=c(1, 2, 1, 1) ),
        main="lh prediction",
        xlab="Time", ylab="lh" )

### fit for deaths data
myplot(deaths)

## stl decomposition of deaths
deaths.stl <- stl(deaths, "periodic")
plot(deaths.stl)

# get the remainder
#rem <- deaths.stl$time.series[,"remainder"]
#trend <- deaths.stl$time.series[,"trend"]
deaths.seasonal <- deaths.stl$time.series[,"seasonal"]

X <- deaths - deaths.seasonal
myplot(X) 
spectrum(X, span=c(3, 3))

deaths.ar1 <- arima(X, order=c(1, 0, 0))
tsdiag(deaths.ar1, gof.lag=30)
deaths.autofit <- auto.arima(X, max.P = 0, max.Q = 0)
tsdiag(deaths.autofit, gof.lag=30)
deaths.autofit <- auto.arima(X)
tsdiag(deaths.autofit, gof.lag=30)


deaths.fit <- auto.arima(X)
tsdiag(deaths.fit, gof.lag=30)


# X  <- rem
deaths.fit <- auto.arima(X)
tsdiag(deaths.fit)
# sur cet exemple il faut mieux laisser la tendance, car on arrive mieux
# a fiter avec


## seasonal ARIMA
deaths.diff <- diff(deaths, 12)
myplot(deaths.diff, lag.max=30)

dd.ar <- ar(deaths.diff)
dd.ar$aic

## ce qui 
deaths.fit1 <- arima(deaths, order=c(2, 0, 0), seasonal=list(order = c(0, 1, 0), 
                 period=12 ) ) 
tsdiag(deaths.fit1, gof.lag=30)


deaths.fit2 <- arima(deaths, order=c(2, 0, 0), seasonal=list(order = c(1, 1, 0), 
                 period=12 ) ) 
tsdiag(deaths.fit2, gof.lag=30)


deaths.fit3 <- arima(deaths, order=c(3, 0, 0), seasonal=list(order = c(1, 1, 0), 
                 period=12 ) ) 
tsdiag(deaths.fit3, gof.lag=30)

# 
deaths.autofit <- auto.arima(deaths)
tsdiag(deaths.autofit, gof.lag=30)


### Nottingam data
# on supprime les deux dernieres annees, pour tester la prediction
nott <- window(nottem, end=c(1936, 12))
myplot(nott)
cat("Il a fait", (min(nott) - 32) / 1.8, "°C en moyenne en Fev 1929 !")
# showlags(nott)

boxplot(split(nott, cycle(nott)), names=month.abb)

nott[110] <- 35 # remove the outlier
nott.stl <- stl(nott, "periodic")
plot(nott.stl)

# on enleve la saison
nott.seasonal <- nott.stl$time.series[, "seasonal"]

nott.stlrem <- nott - nott.seasonal
# nott.rem <- nott.stl$time.series[, "remainder"]
# nott.rem <- nott.stlrem - mean(nott.stlrem)
myplot(nott.stlrem)
ar(nott.stlrem)$aic

## on fitte avec un modele AR(1)
nott.ar1 <- arima(nott.stlrem, order=c(1, 0, 0))
quartz()
tsdiag(nott.ar1, gof.lag=30)

nott.fore1 <- predict(nott.ar1, 36)
# mais il faut remettre la saison et la moyenne
#nott.fore1$pred <- nott.fore1$pred + mean(nott.stlrem) + nott.seas[1:36]

nott.fore1$pred <- nott.fore1$pred + nott.seasonal[1:36]

quartz()
ts.plot(window(nottem, 1937), nott.fore1$pred, 
           nott.fore1$pred - 1.96*nott.fore1$se, nott.fore1$pred + 1.96*nott.fore1$se, 
            gpars=list( lty=c(1, 1, 2, 2), col=c("black", "red", "blue", "blue"), 
             lwd=c(2, 2, 1, 1), ylim=c(35, 73) ) )
legend( "topleft", c("True value", "Predition", "95% Confidence"), 
   lty=c(1, 1, 2, 2), lwd=c(2, 2, 1, 1),  col=c("black", "red", "blue", "blue") )

# Autre approche en utilisant la differentiation
nott.diff <- diff(nott, 12)
myplot(nott.diff, lag.max=30)
spectrum(nott.diff)

nott.fit1 <- arima(nott, order=c(1, 0, 0), seasonal=list( order=c(2, 1, 0), period=12 ) )
nott.fit1
tsdiag(nott.fit1, gof.lag=30)

nott.fit2 <- arima(nott, order=c(0, 0, 2), seasonal=list( order=c(0, 1, 2), period=12 ) )
nott.fit2
tsdiag(nott.fit2, gof.lag=30)

nott.fit3 <- arima(nott, order=c(1, 0, 0), seasonal=list( order=c(0, 1, 2), period=12 ) )
nott.fit3
tsdiag(nott.fit3, gof.lag=30)

nott.fit4 <- arima(nott, order=c(1, 0, 2), seasonal=list( order=c(2, 1, 2), period=12 ) )
nott.fit4
tsdiag(nott.fit4, gof.lag=30)

nott.autofit <- auto.arima(nott)
nott.autofit
tsdiag(nott.autofit, gof.lag=30)

nott.fore2 <- predict(nott.fit4, 36)


# comparaison entre l'approche STL et l'approche SARIMA

layout(1:2)

ts.plot(window(nottem, 1937), nott.fore1$pred, 
           nott.fore1$pred - 1.96*nott.fore1$se, nott.fore1$pred + 1.96*nott.fore1$se, 
            gpars=list( lty=c(1, 1, 2, 2), col=c("black", "red", "blue", "blue"), 
             lwd=c(2, 2, 1, 1), ylim=c(35, 73) ) )
title("Using STL")

legend( "topleft", c("True value", "Predition", "95% Confidence"), 
   lty=c(1, 1, 2, 2), lwd=c(2, 2, 1, 1),  col=c("black", "red", "blue", "blue") )

ts.plot(window(nottem, 1937), nott.fore2$pred,
           nott.fore2$pred - 1.96*nott.fore1$se, nott.fore2$pred + 1.96*nott.fore2$se, 
            gpars=list( lty=c(1, 1, 2, 2), col=c("black", "red", "blue", "blue"), 
             lwd=c(2, 2, 1, 1), ylim=c(35, 73) ) )
title("Using Box-Jenkins")
legend( "topleft", c("True value", "Predition", "95% Confidence"), 
   lty=c(1, 1, 2, 2), lwd=c(2, 2, 1, 1),  col=c("black", "red", "blue", "blue") )


# on peut refaire la meme chose avec
nott.ar1 <- auto.arima(nott.rem)
# mais cela n'ameliore pas tellement la prediction




## Donnes AirPassenger
ap <- AirPassengers
log.ap <- log(ap)
log.ap.stl <- stl(log.ap, "periodic")
log.ap.trend <- log.ap.stl$time.series[, "trend"]
X <- diff(log.ap - log.ap.trend, lag=12)
myplot(X, lag.max=30)




arima(X, order = c(0, 0, 1), seasonal=list( order=c(1, 0, 1), period=12) )

arima(log.ap - log.ap.trend, order = c(0, 0, 1), seasonal=list( order=c(1, 0, 1), period=12))


log.ap.autofit <- auto.arima(log.ap)
tsdiag(log.ap.autofit, gof.lag=30)



log.ap.fit1 <- arima(X, order = c(1, 12, 1), seasonal=list( order=c(1, 0, 1), period=12) )
tsdiag(log.ap.fit1, gof.lag=30)

X.autofit <- auto.arima(X)
tsdiag(X.autofit, gof.lag=30)

ap.log <- log(ap)
myplot(ap.log)


